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1  Introduction 


Scintillation  is  a  disturbance  in  the  amplitude  and  phase  of  a  radio  wave  caused  by  random 
electron  densities  in  the  ionosphere  [1].  The  electron  densities  along  the  path  of  a  propagating 
radio  signal  cause  the  scintillation  [2],  These  electron  densities  are  a  result  of  disturbances  in  the 
ionosphere  and  because  the  densities  are  the  driving  force  behind  scintillation,  the  ensuing 
change  in  the  signal’s  phase  is  random  as  well  [3].  Scintillation,  especially  around  the  equator, 
interferes  with  Air  Force  satellite  communication  systems  as  well  as  the  Global  Positioning 
System  (GPS).  Disturbances  in  the  ionosphere  are  dependent  on  the  season,  solar  cycle,  and 
other  variable  factors  [4].  Depending  on  the  severity  of  the  ionospheric  disturbance,  scintillation 
will  be  more  or  less  pronounced. 

The  ionospheric  scintillation  problem  can  be  considered  as  a  signal  propagation  problem 
with  the  perturbations  in  ionosphere  depicted  as  one  or  multiple  random  phase  screens.  These 
phase  screens  cause  small  localized  perturbations  in  the  signal  phase,  which  will  cause  amplitude 
perturbations  due  to  diffraction  and  localized  interference.  A  popular  computational  technique  for 
modeling  ionospheric  scintillation  is  the  phase  screen  method.  It  calculates  monochromatic  plane 
wave  propagation  through  any  given  phase  screen  pattern  and  predicts  amplitude  and  phase 
perturbations.  Signals  that  are  not  plane  waves  are  analyzed  by  Fourier  transforming  the  signal 
and  perfonning  calculations  at  component  frequencies.  Although  computationally  efficient,  the 
phase  screen  method  is  limited  to  simulating  one  frequency  at  a  time  and  is  confined  to  relatively 
small  angles  [5].  In  other  words,  this  method  is  used  to  model  a  plane  wave  that  does  not  deviate 
substantially  from  its  original  direction  of  propagation. 

Recent  advancements  in  computational  capabilities  provide  new  means  of  modeling  the 
wave  propagating  through  the  ionosphere  and  permit  analysis  of  scintillation  characteristics  in  a 
more  detailed  manner.  One  such  method  is  the  Finite  Difference  Time  Domain  (FDTD)  method, 
pioneered  by  Kane  Yee  in  1966  [6].  FDTD  is  a  grid  based,  full-wave  technique  used  to  solve 
Maxwell’s  differential  equations  in  the  time  domain.  Kane  Yee  developed  a  set  of  finite 
difference  equations  that  are  central-difference  in  both  time  and  space  and  simultaneously 
encompass  Maxwell’s  equations  on  a  microscopic  and  macroscopic  level  [7].  Yee’s  algorithm 
provides  a  smooth  transition  to  a  discretized  fonn  of  the  relationship  between  the  electric  (E)  and 
magnetic  (H)  fields  by  utilizing  the  coupled  nature  of  E  and  H  given  in  Maxwell’s  equations. 
The  robustness  of  this  technique  stems  from  using  the  information  from  both  the  E  and  H  fields. 
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By  using  both  fields,  different  modes,  material  properties,  and  features  unique  to  each  field  are 
modeled  through  straightforward  variations  in  the  FDTD  algorithm  [7], 

The  goal  of  this  research  activity  was  to  develop  2-D  FDTD  simulations  for  a  better 
representation  of  the  ionospheric  scintillation  effects  in  communications.  FDTD  is  versatile 
enough  to  enable  numerical  experiments  on  various  parts  of  the  model  (ionosphere,  transmission 
sources,  communication  signal  specifics,  etc.)  and  can  help  in  fonnulation  of  solutions  that 
minimize  scintillation  effects. 

Section  2  of  this  report  provides  an  overview  of  the  numerical  methods.  Section  3 
discusses  issues  specific  to  FDTD  implementation  of  the  ionospheric  scintillation  problem. 
Section  4  presents  results  obtained  and  lessons  learned.  Section  5  provides  recommendations  for 
future  research  directions.  Section  6  concludes  the  report. 
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2  Modeling  and  Methods 

2.1  Phase  Screen  Method 

The  phase  screen  method  models  the  ionosphere  as  a  change  in  the  incident  wave’s  phase 
along  each  point  of  the  wave  front  [1].  The  number  of  points  depend  on  the  resolution  of  the 
chosen  model.  For  example,  the  sample  1-D  phase  screen  output"  provided  by  the  Air  Force 
Research  Laboratory  (AFRL)  for  this  research  is  a  20  kilometer  wide  distribution  composed  of 
20,000  points.  Therefore,  the  screen  contains  one  point  per  meter,  which  results  in  a  sampling 
resolution  of  1  meter.  The  electron  densities  in  the  ionosphere  fonn  striated  patterns  because  of 
their  tendency  to  line  up  with  earth’s  magnetic  field  [3].  This  alignment  is  important  because  it 
leads  to  a  dimensional  simplification  in  the  ionospheric  modeling.  Near  the  equator  a  1-D  phase 
screen  is  accurate  because  of  the  earth’s  horizontal  magnetic  field  and  because  the  ionosphere 
has  a  relatively  insignificant  thickness  [1]. 

The  phase  screen  method  [8-10]  is  based  on  the  parabolic  wave  equation.  This  is  an 
approximation  of  the  wave  equation  where  the  solution  is  weakly  perturbed  to  deviate  from  a 
straight  propagation  direction.  Starting  from  the  2-D  Helmholtz  equation 


d2xp  d2xp 


(1) 


a  7  +  TT  +  k  n  Ip  =  0 
ox1  dz z 

and  assuming  slowly  changing  refractive  index  n,  a  harmonic  time  dependence  e-twt,  no 
changes  in  the  y  direction,  and  a  perturbed  solution 

u(x,z )  =  e~lkxip(x,z)  ^ 2 ^ 

one  can  insert  the  prescribed  solution  into  the  equation  to  get 

d2u  d2u  du  ,  . 

— ■  +  +  2i/c— +  /c2(n2  -  l)u  =  0  (3) 

ox1  dz 1  dx 


This  equation  can  be  factorized  as 

(Ia  ik(i  -  Q'>}{^+ik(i + «}“ = ° 

2  Pedersen,  T.,  private  communication 
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with  Q  —  J^2^2+n'2  an  operator.  Equating  the  first  term  to  zero,  one  gets 

g+ik(1.Q)Ju  =  0  (5, 

Assuming  a  refractive  index  close  to  air  (n  ~  1)  and  a  plane-wave  like  solution  propagating  at 
small  angles  near  the  x  direction,  one  can  approximate  Q  and  the  square  root  in  its  definition  to 
end  up  with  a  parabolic  differential  equation  which  is  second  degree  in  z  and  first  degree  in  x, 


du 

2  ik—  = 
ox 


dz2 


+  k2(n2  —  1)  [it 


(6) 


The  phase  screen  method  identifies  irregularities  in  the  ionosphere  as  layers  of  phase  screens. 
Propagation  through  a  phase  screen  piece  or  a  section  of  free  space  is  then  calculated  via 
analytical  solutions  of  the  parabolic  differential  equation. 

Advantages  of  this  method  include  high  efficiency  in  calculating  scintillation  effects. 
Disadvantages  are  due  to  limitations  in  constructing  the  parabolic  differential  equation,  such  as 
being  limited  in  bandwidth  (constant  or  near-constant  frequency)  and  propagation  angle, 
assumption  of  plane  wave-like  solutions  and  oversimplified  representation  of  ionospheric 
disturbances. 


2.2  Finite  Difference  Time  Domain  (FDTD)  Method 

The  FDTD  method  is  based  on  a  rather  direct  approach,  where  partial  derivatives  in  space  and 
time  in  half  of  Maxwell’s  equations  are  expressed  as  second-order  accurate  finite  differences  on 
a  staggered  grid  of  electric  and  magnetic  fields.  Maxwell’s  curl  equations  (equations  1  and  2)  are 
split  into  x,  y,  and  z  components.  The  result  is  six  coupled  scalar  equations  for  six  field 
components. 


dH  _ 
dt 

dE 

dt 


1  1 

--V  XE--M 

p.  p 

1 

=  -VxH-J 

£ 


(7) 

(8) 


For  2-D  implementation,  differentiation  with  respect  to  z  is  set  equal  to  zero  and  the  six  coupled 
equations  reduce  to  the  following  scalar  fonn  (from  [7]): 
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djh 

dt 

dH„ 


dJk 

fi  L  dy 
dEv 


dt 


li  L  dx 


dHz  _ 

1 

r dEx 

\ 

1 

dt 

E 

.  dy 

dx 

-  (Mx  +  a*Hx ) 

(My  +  0*Hy) 


-  (Mz  +  a*Hz ) 


dE, 

dt 

dEy 

dt 


1  \dH7 


sidy 
dH7 


dx 


dEz 

1 

r  dHy 

dt 

£ 

.  dx 

dy 

( Jx  "I"  oEx) 

(Jy  +  OEy) 

( /z  +  0-^z) 


(9) 

(10) 

(11) 

(12) 

(13) 

(14) 


where  x,  y  and  z  are  spatial  coordinate  axes,  t  denotes  time,  s  and  y  denote  pennittivity  and 
penneability,  a  is  electric  conductivity,  cr*  is  magnetic  loss,  Hx,  Hy  and  Hz  are  components  of 
magnetic  field  strength,  Ex,  Ey  and  Ez  are  electric  field  components,  Jx,  Jv  and  J:  are  current 
density  components,  and  Mx,  My  and  M-  are  components  of  magnetization.  These  equations  are 
the  basis  of  the  FDTD  algorithm  in  2-D.  The  fields  in  2-D  are  uncoupled  into  two  orthogonal 
representations.  Equations  (9),  (10)  and  (14)  constitute  the  TMZ  mode  and  equations  (11),  (12) 
and  (13)  the  TEZ  mode.  A  similar  dimensional  reduction  is  performed  for  a  1-D  implementation, 
which  will  not  be  covered  here. 

If  one  adopts  the  following  convention  for  notational  convenience 


(ij)  =  (iAx.jAy)  15 ) 

u(iAx,jAy,nAt )  =  ufj 

where  Ax,  Ay  are  the  spatial  discretization  units  and  At  is  the  unit  increment  in  time.  The  first 
partial  space  derivative  of  the  entity  u  in  x  is  then  written  as 

du(iAx,jAy,nAt)  Ui+\.i 

dx  Ax 

Similarly,  the  time  derivative  is  written  as 

n+i  n-i 

du(iAx,jAy,  nAt)  _  \j  ~  \j  (17) 

dt  At 

The  discretization  of  the  electric  and  magnetic  field  components  is  perfonned  in  such  a  way  that 
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they  are  staggered  in  space  (Figure  1)  as  well  as  time  (i.e.  if  electric  field  components  are 
expressed  at  time  points  (n  +  ^  At  and  (n  —  At,  magnetic  field  components  are  expressed  at 
nAt  and  (n  +  l)At).  At  this  point,  one  can  choose  to  collect  terms  in  equations  such  that  fields  at 
time  points  (n  +  ^  At  and  (n  +  l)At  are  expressed  in  terms  of  fields  at  earlier  points  in  time 

(e.g.  (n  —  At,  nAt,. . .).  This  results  in  an  explicit  marching  scheme. 


Hr(ie-\,je) 

' 

o 

'  O  f  o 

I  °  I  °  1  ° 

1  ° 

I  ° 

f  E.Me,je) 

f/,(l,ye-l) - ► 

—  — 

t —  t — *  t — * 

♦  * 

1 _ 

♦  * 

_ Hx(ie,je- 1) 

O 

o  |  o 

O  O  [  o 

O 

o 

O 

»  t  * 

K  t  *  t  * 

1 

1 _ 

4  * 

- *. 

o 

o  |  o 

O  O  T  o 

fT 

o 

o 

f  *  t  * 

w  • _  l 

K  t  *  t  * 

1 

A  *  A 

t — ft. 

o 

o  f  o 

O  o  1  o 

I"®" 

O 

o 

1  _  1 

►  %  *■- 

K  *  t  *  ♦  * 

t  — - ♦ 

1  _  1 

f  *4 

o 

o  }  o 

0 

0 

0 

O 

o 

o 

fi8(U) 

1  1 
( 

• 

• 

t 

1  1 

1 

1  , 

j  i 

Ez(ie,  1) 

£t(2,l) 

Ez{ie- 

1.1)  1 

Figure  1  2-D  FDTD  grid  with  TMZ  components1 


Derivation  of  Yee’s  algorithm  will  not  be  covered  here,  as  it  can  be  found  in  many 
resources,  for  example  see  [7],  The  result  for  Ez,  after  discretization  of  equation  (14)  is 
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Since  the  electric  and  magnetic  field  terms  are  staggered  in  time  and  expressed  in  tenns  of  one 
another,  the  explicit  marching  scheme  is  invoked  by  calculating  electric  and  magnetic  field 
components  one  after  another,  while  the  time,  nAt,  is  increased  until  a  predetermined  value  or  a 
simulation  state  is  reached.  Spatial  staggering  requires  calculation  of  field  components  on  the 
entire  grid  at  every  time  step,  as  the  calculation  formulas  refer  to  field  values  at  neighboring  grid 
locations.  This  also  means  field  values  at  grid  boundaries  cannot  be  calculated  with  the  usual 
formulas  and  must  be  determined  in  other  ways  (i.e.  by  formulating  boundary  conditions). 

Entities  referring  to  material  properties  such  as  pennittivity,  permeability  and 
conductivity  are  discretized  appropriately  and  utilized  in  defining  the  geometry  of  objects  on  the 
grid.  Other  entities  such  as  ad-hoc  source  terms  J  and  M  require  special  attention  as  they  provide 
both  means  of  convenience  (e.g.  one  can  inject  energy  into  the  FDTD  simulation  by  providing  a 
current  J  at  an  appropriate  location)  and  trouble  (e.g.  injected  currents  are  unphysical  and  such 
energy  sources  will  cause  spurious  reflections  when  certain  conditions  are  met). 

One  particular  source  type  used  throughout  this  project  is  known  as  Total  Field/Scattered 
Field  (TF/SF)  source.  The  TF/SF  technique  is  advantageous  because  it  implements  plane  wave 
propagation  with  minimal  spurious  reflections.  The  foundation  of  the  TF/SF  technique  is  its 
dependence  on  the  linearity  of  Maxwell’s  equations.  The  FDTD  simulation  is  divided  into  two 
regions  on  a  closed  virtual  boundary,  named  total  field  and  scattered  field  regions.  The  total  field 
region  is  enclosed  by  the  scattered  field  region  (Figure  2).  The  following  relationships  are 
asserted  on  the  virtual  boundary  separating  the  two  regions 


E  total  Einc  +  E  scat- 
H total  —  El inc  T  H scat 


(19) 

(20) 


where  Einc  and  Hinc  represent  the  values  of  the  incident  fields  and  Escat  and  Hscat  represent 
the  values  of  the  scattered  fields.  The  incident  field  is  introduced  through  the  virtual  boundary 
via  special  field  updates  such  that  incident  fields  are  present  in  the  total  field  region  but  not  in  the 
scattered  field  region. 
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Truncation 

Figure  2  TF/SF  space  lattice  depicting  the  total  field  and  scattered  field  regions 


The  special  updates  required  for  implementation  of  the  TF/SF  technique  are  achieved  by 
fixing  inconsistencies  between  regions.  The  updates  are  best  described  with  a  1-D  example  as 
shown  in  Figure  3.  In  the  1-D  case  the  virtual  boundary  between  regions  becomes  a  single  point 
on  each  side  of  the  total  field  region. 
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Figure  3  1-D  TF/SF  space  lattice  showing  components  in  need  of  special  updates1 


At  the  left  interface  between  regions  1  and  2,  Ez  totai  will  have  an  inconsistent  update  [7] 

^  1 71+1  ^  1 71  i  At  |n+J  IT  |n+j\  ^ 

Ez, total  |;  —  ^z, total  L-  T  —  I  Hy  total  \ .  1“  ^y,scat|.  l)  (71) 

L  L  SqlaX  y  '  2  2/ 

where  iL  is  the  coordinate  of  the  left  TF/SF  boundary.  Clearly,  a  total  field  component  cannot  be 
updated  by  a  scattered  field  component.  To  correct,  the  incident  field  term  must  be  added  to  the 
scattered  field  term 
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^  |J1+1  ^  |  n  1  At  |n+J  rr  |  n  +  J  TJ  |n  +  j'\  /  TT  x 

Ez,total\ :  —  Ez  total]:  T  ~  T~  l  Hy  total  \ .  1  ^v,scat  .  l~^y,inc\.  1  I  ( -- ) 

L  L  SqLAJC  y  I'll'  2  2  ^  2/ 

where  Eztotal  is  now  updated  only  by  total  fields.  A  similar  inconsistency  around  the  right 
boundary  at  iR  exists 
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and  it  is  handled  in  a  similar  way 


.n+l  .n  At  /  |n+2  \n+2  \n+2  \ 

Ez, total  |;  —  Ez  total  |;  ^  7  I  Ely  scat  | .  l  —  ^v, total  | .  ld”^v,tnc|.  1  I  (24) 

LR  LR  EqIaX  \  lR  +  2  1r~2  1r  +  2 ) 


but  one  must  note  the  swapped  position  of  the  total  and  scattered  field  terms  and  the  resulting 
sign  reversal  of  the  incident  field  tenn. 

A  numerical  method  is  considered  numerically  stable  if  it  does  not  have  any  numerical 
artifacts  to  cause  monotonic  growing  or  shrinking  of  the  calculated  tenns.  For  the  explicit 
marching  FDTD  formulation,  simulation  stability  is  dependent  on  the  size  of  the  time  step. 
Detailed  analysis,  found  elsewhere  [7,11],  leads  to  a  relationship  between  the  time  step  size  (At) 
and  the  grid  cell  size  (Ax,  Ay)  known  as  the  Courant  stability  limit: 


1 

At  <  - 
c 

where  c  is  the  speed  of  light  or  electromagnetic  energy  propagation. 

A  related  concept,  numerical  accuracy,  refers  to  the  finite  accuracy  of  modeling  in 
numerical  simulations.  For  FDTD,  it  means  the  model  must  be  able  to  adequately  resolve  both 
the  smallest  physical  feature  and  the  smallest  wavelength  in  the  simulation.  While  the  Nyquist- 
Shannon  sampling  theorem  [12]  requires  at  least  two  mesh  cells  per  physical  feature  or 
wavelength,  common  practice  lore  favors  usage  of  better  resolution,  where  8-10  mesh  cells  is 
considered  a  good  number.  It  must  be  noted  that  increasing  resolution  in  an  explicit  marching 
based  FDTD  scheme  will  require  reduction  of  the  time  step  size  per  the  Courant  stability  limit. 
FDTD  offers  many  advantages  for  ionospheric  scintillation  modeling,  such  as 


1  +  1 


I  (Ax)2  1  (Ay)2 


(25) 
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•  Having  a  full-wave  simulation  not  limited  by  the  parabolic  differential  equation 
or  any  approximations  in  bandwidth,  propagation  angle  and  source  wave  forms 

•  The  possibility  of  formulating  dispersive,  absorptive,  gain,  plasma  or  other  types 
of  media  for  better  representation  of  the  ionosphere  and  sources  of  scintillation 

•  The  ability  to  investigate  transients  and  time-dependent  phenomena  (unavailable 
when  a  frequency-domain  method  is  used) 

•  The  ability  to  perform  rigorous  numerical  experiments  beyond  insight-gathering 
phase 

•  High-performance  implementation  possibilities  to  keep  simulation  times  at  a 
manageable  level 

The  disadvantages  of  FDTD  include  increased  requirements  for  computing  resources  and  time, 
intricacy  of  implementation,  numerical  dispersion  and  other  numerical  artifacts  for  large  and/or 
oblong  grid  configurations,  difficulty  modeling  curved  shapes  without  approximations  and 
difficulty  using  grids  with  many  different  element  sizes  and/or  shapes. 

2.3  Near-to-Near  and  Near-to-Far  Field  Transforms 

Near-to-near  field  and  near-to-far  field  transforms  were  developed  to  address  a  crucial  deficiency 
in  FDTD.  Due  to  the  way  FDTD  is  fonnulated,  wave  propagation  speed  on  an  FDTD  mesh 
changes  slightly  depending  on  propagation  direction  (a  phenomenon  known  as  numerical 
dispersion  [7]).  This  means  a  wave  front  propagating  in  multiple  directions  in  a  large  FDTD  grid 
becomes  increasingly  distorted  as  it  continues  to  propagate.  Near-to-near  field  and  near-to-far 
field  transforms  use  different  techniques  to  extrapolate  electromagnetic  field  values  to  locations 
of  interest  without  having  to  extend  the  entire  FDTD  grid  to  cover  such  locations. 

One  commonality  among  all  such  transforms  is  their  dependence  on  the  surface 
equivalence  theorem,  which  itself  is  a  reformulation  of  Huygens’  principle.  Surface  equivalence 
theorem  establishes  the  equivalence  of  electromagnetic  fields  generated  by  sources  and 
scattering  inside  a  closed  surface  with  electromagnetic  fields  generated  by  appropriate  electric 
and  magnetic  current  distributions  on  the  closed  surface.  The  surface  does  not  have  to  be  a 
physical  object  and  can  be  in  any  convenient  shape.  In  the  case  of  near-to-near  field  and  near-to- 
far  field  transforms,  surface  equivalence  theorem  serves  to  replace  electromagnetic  fields 
generated  inside  a  virtual  boundary  (Figure  4(a))  with  equivalent  surface  currents  on  the  surface 
(Figure  4(b)).  This  simplifies  calculation  of  a  number  of  entities  such  as  the  far  field  response, 
since  the  inside  of  the  virtual  boundary  is  devoid  of  electromagnetic  fields  and  any  objects  inside 
it  can  be  modified  or  removed  without  affecting  the  calculation  results  of  interest  (Figure  4(c)).  If 
free  space  is  present  outside  the  virtual  boundary,  it  becomes  convenient  to  assume  free  space  is 
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present  inside  the  boundary.  Then  any  calculations  of  interest  can  be  assumed  to  take  place  in 
free  space,  devoid  of  scatterers. 


Figure  4  (a)  A  virtual  closed  boundary  around  a  scatterer  or  source,  shown  with  fields  across  the  boundary,  (b)  Fields  inside  the 
boundary  are  set  to  zero.  Fields  outside  the  boundary  are  generated  by  surface  currents  on  the  boundary,  (c)  The  scatterer  or 
source  can  now  be  disposed  of  since  its  presence  or  absence  is  irrelevant  in  a  region  devoid  of  electromagnetic  fields.  From  [13]. 


There  are  a  number  of  near-to-near  field  and  near-to-far  field  transformation  methods 
formulated  for  FDTD  in  literature,  based  on  integration  of  surface  currents  [7,14-22],  multipole 
expansion  [23,24],  Fresnel  diffraction  integral  [25]  and  Kirchhoff’s  surface  integral  [26].  Of 
these  methods,  two  were  employed  in  the  course  of  this  research:  The  time  domain  near-to-far 
field  transfonn  by  Luebbers  et  al.  [14]  and  Kirchhoff’s  surface  integral  representation  (KSIR) 
based  near-to-near  field  transform  by  Ramahi  [26].  These  methods  will  be  discussed  below. 
Further  implementation  details  specific  to  the  ionospheric  scintillation  problem  can  be  found  in 
Section  3. 

The  most  fundamental  equations  to  the  time  domain  near-to-far  field  transfonn  by 
Luebbers  et  al.  [14]  are  given  by 


W(r,t)  = 

U(r,t)  = 
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where  Eg  and  E(f)  are  the  spherical  far  field  E  field  components,  W  and  U  are  the  far-field 
potential  vectors,  Js  and  Ms  are  the  electric  and  magnetic  currents  along  the  virtual  surface,  S, 

r — r'-r 

and  are  expressed  in  retarded  time,  t - - — .  This  term  represents  the  time  delay  between  the  E 

and  H  fields  on  S  to  their  appearance  at  the  far  field  point.  The  surface  S  should  cover  all  sources 
in  the  FDTD  simulation.  If  a  total  field/scattered  field  source  is  used,  S  will  be  in  the  scattered 
field  (i.e.  outer)  region.  A  typical  configuration  for  a  2-D  FDTD  simulation  is  depicted  in  Figure 
5. 


Figure  5  Typical  configuration  of  a  2-D  FDTD  simulation  utilizing  near-to-far  field  (NTFF)  transformation  and  total 
field/scattered  field  (TFSF)  source.  PML  refers  to  a  perfectly  matched  layer  absorbing  medium  placed  on  the  outer  sides  of  a 
scattering  simulation  to  prevent  reflection  of  outgoing  energy. 

The  general  procedure  for  the  near-to-far  field  transformation  is  to  evaluate  the  integrals 
in  equations  (26)  and  (27).  To  give  an  example,  for  the  surface  excitation  Mz,  located  within  a 
rectangular  patch  AxAz,  the  contribution  from  the  integral  in  equation  (27)  can  be  written  as 


1  d  AxAz  d  ,  s 

AU  —  AUzz  —  - - —  (MzzAxAz)  =  - - —  (Exz)  < 30 ) 

4nrc dt  4nrc  dt v  J 

The  number  of  time  steps  until  AU  from  contributes  to  the  far  field  vector  potentials  is  given  by 

=  r-r  -f  (31) 

J  cAt 

which  is  the  total  time  delay  divided  by  the  FDTD  time  step.  Thus,  the  integral  in  equation  (27) 
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is  written  in  finite  difference  notation  as 


n+^+f  ,n+h+f 

Uz\r  2  =  Uz\r  2  + 


Ax  A z  d 
4nrc  dt 


(Ex  I 


n+ 1 
xi  n 


Ex\Ur 0 


(32) 


1 

It  is  highly  unlikely  that  n  +  -  +  f  will  be  an  integer.  This  means  the  near-to-far  field 

transformation  will  be  incompatible  with  the  general  FDTD  semantics.  A  solution  is  to  make  two 
fractional  adjustments  based  on  linear  interpolation  to  the  nearest  integer  time  steps,  given  as 


AxAz  d 

uz\r  =  uzr  +  a  - 

„  ,  „  AxAz  d 

tt  inn+1  _  it  inn+1  i  „  .  rr?  in+1  rr  in'v 

uz\r  ~  uz  Ir  +  a  ^nrc  n  ~  Lx\  n) 

1 

where  nn  is  the  integer  truncation  of  the  value  nt - 1 -  f  and  a  is  given  by 

a  =  (n+\  +  f)-nn 


(33) 

(34) 


(35) 


which  represents  the  fractional  time  step  between  the  exact  delay  of  A  U  and  the  integer  time  step 
that  is  just  prior.  Thus,  (1  —  a)  is  representative  of  the  fractional  time  step  between  the  exact 
delay  of  A  U  and  the  integer  time  step  that  is  just  after  it. 

The  transformation  method  outlined  above  is  valid  only  if  the  location  of  interest  is  in  the 
far  field  region  of  the  aperture  under  study.  This  requirement  is  not  always  satisfied  for  the 
ionospheric  scintillation  problem,  especially  if  one  wishes  to  compare  FDTD  results  against 
phase  screen  method  results  due  to  the  size  of  the  phase  screen  and  the  distance  from  it.  As  such, 
a  more  general  method  that  can  handle  both  near  and  far  field  extrapolation  is  needed. 
Kirchhoff ’s  surface  integral  representation  (KSIR)  based  near-to-near  field  transfonn  by  Ramahi 
[26]  is  one  such  method  and  was  implemented  during  the  course  of  this  research  project. 

KSIR  is  based  on  Kirchhoflf’s  surface  integral,  which  relates  a  field  y/  inside  a  closed 
volume  V  to  the  field  and  its  derivatives  on  the  surface  of  V, 


V'ip(x',t')_R_  ,  ,  R  dip(xr ,  t') 

r>3  f  z 


da' 


(36) 


Jret 


R  l?3  T  cR2  dt ' 

where  R  —  x  —  x',  R  =  \r\,  n  is  the  unit  normal  vector  to  the  surface,  c  is  the  speed  of  light  and 
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the  subscript  ret  indicates  the  integral  is  evaluated  at  the  retarded  time  t'  =  t  —  R/c.  Ramahi 
goes  on  to  provide  part  of  a  3-D  FDTD  derivation,  where  xp  —  Ex  evaluated  at  FDTD  coordinate 
index  k  —  k0  of  an  xv-plane  surface,  t*  —  t  +  R/c  is  the  advanced  time  (in  lieu  of  the  retarded 
time  in  the  earlier  fonnula),  and  a  contribution  to  the  integral  Ex  k  (x,  t/J  is  given  by 


Ex,k0(x>  tn)  =  Fiin)  +  F2(n  +  1)  +  F3(n  +  2) 


dE* 

dz 


Fi(n)  =  ^^;Ex(i',j',k0,n)  Aj 


F2(n  +  1)  =  ADZ  +  B)Ex(i',j',k0,n  +  l)Ai/j/ 


ZC 

—  Ex(i',j',k0,n  +  2)Ai/j/ 
i'i' 

1  1 

A  = - 

47T  R 

1  —  cos  0' 


47T  R3 
1  —cos  6' 


An  cR2 

=  DzEx(i,j,  k0,n  +  1) 

_  Ex(i,j,k0  +  l,n  +  1)  -  Ex(i,j,k0  -  l,n  +  1) 
“  2Az 


(37) 

(38) 

(39) 

(40) 

(41) 

(42) 

(43) 

(44) 


where  the  primed  indices  belong  to  the  KSIR  integration  (vs.  the  unprimed  indices  referring  to 
FDTD),  Aj/j/  is  the  area  of  a  subsurface  over  which  the  integration  is  perfonned,  R  is  the  distance 
and  6'  is  the  angle  that  the  normal  of  subsurface  Aj/p  makes  with  the  point  of  interest  x.  The 
interpretation  of  equation  (37)  is  as  follows:  At  FDTD  time  step  n,  the  terms  F1(n),  F2(n)  and 
F3(n)  are  computed,  and  the  contribution  F{  (n)  is  added  to  Ex(n*),  where  n*  —  |(n+  1)  + 
R/(cAt)\int  is  the  propagation  delay  in  terms  of  FDTD  time  steps,  truncated  to  the  nearest 
integer  not  exceeding  its  value.  Similarly,  at  FDTD  time  step  n  +  1,  the  tenn  F2(n  +  1) 
contributes  to  Ex(n *),  and  at  FDTD  time  step  n  +  2,  the  term  F3(n  +  2)  contributes  to  the  same 
term.  Thus,  at  every  FDTD  time  step,  three  contributions  per  subsurface  App  are  made  to 
different  Ex  terms  with  different  propagation  delays. 
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3  FDTD  Implementation  of  the  Ionospheric  Scintillation 
Problem 

3.1  Ionosphere  as  a  Phase  Screen 

One  of  the  goals  of  this  research  project  was  to  create  a  2-D  FDTD  simulation  of  a  1-D 
phase  screen  and  compare  the  results  to  a  phase  screen  method  based  simulation.  The  phase 
screen  simulation  is  representative  of  a  plane  wave  traveling  through  the  equatorial  ionosphere 
and  impinging  Earth’s  surface.  Thus  the  goal  of  the  FDTD  simulation  was  to  replicate  the  phase, 
amplitude,  and  amplitude  scintillation  index  (S4)  provided  by  the  phase  screen  simulation  at  both 
the  ionosphere  and  earth’s  surface.  The  amplitude  scintillation  index  is  defined  in  terms  of  root- 
mean-square  variance  of  the  irradiance  /  due  to  a  wave  [1] 

S%  —  «/2)  —  (i)2)/(i)2  (45) 

where  the  brackets  indicate  spatial  or  temporal  average.  In  order  to  replicate  the  phase  screen 
method  scenario  in  FDTD,  a  large  grid  with  high  resolution  (i.e.  small  grid  cells)  and  an  incident 
plane  wave  are  required. 

The  fluctuations  in  ionosphere  were  implemented  in  FDTD  as  a  dielectric  slab  that 
represents  the  ionosphere  as  a  simple  dielectric  with  varying  thicknesses  intended  to  create 
predetennined  phase  delays.  A  plane  wave  source  injected  over  a  total  field/scattered  field  source 
boundary  propagated  through  the  dielectric.  The  electric  field  phase  was  measured  and  compared 
to  a  reference  case  with  a  constant-thickness  dielectric  slab  or  no  slab  to  obtain  phase  delay 
values.  The  phase  delay  values  were  calculated  with  and  without  a  near-to-near  field 
transformation  based  on  Kirchhofif’s  surface  integral  representation  (KSIR).  In  the  latter,  the 
FDTD  simulation  grid  was  elongated  slightly  to  cover  the  line  of  extrapolation  provided  by  the 
near-to-near  field  transfonnation  in  the  fonner  case. 

The  phase  screen  data  provided  by  AFRL1  for  this  project  was  a  20  km- wide  distribution 
of  ionosphere  with  a  sampling  resolution  of  1  m.  The  phase  screen  contained  both  positive  and 
negative  phase  values.  Positive  phase  values  cannot  be  represented  by  a  length  of  dielectric  slab 
because  a  wave  propagating  through  a  dielectric  only  experiences  phase  delay  in  comparison  to 
free  space.  Therefore,  every  phase  point  is  referenced  to  the  greatest  of  the  phase  points;  the 
greatest  phase  value  would  then  have  zero  phase  delay.  The  new  reference  is  implemented 
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according  to  the  following  equation 


* Prefix. )  —  (p{x )  (Pmax 

where  (Pmax  is  the  largest  phase  value,  (()  is  the  original  array  of  phase  screen  values,  and  (pref  is 
the  newly  referenced  phase  screen  array.  The  variables  are  a  function  of  x,  which  represents 
position  along  the  ionosphere.  By  applying  this  equation,  the  maximum  delay  is  now  equal  to 
zero  and  every  other  phase  is  certainly  negative,  i.e.  delayed.  This  procedure  is  merely  the 
equivalent  of  shifting  the  original  phase  screen  phase  values  down  by  a  constant  (pmax-  After  the 
phases  have  been  referenced,  the  slab  thickness  for  every  phase  is  calculated  as  follows: 


](r}  —  Prefix') 

[-27 zfiJe^Tr  -  1  XV^oMo)]  (47) 

where  l(x )  is  the  length  corresponding  to  a  specific  phase  delay,  /  =  250  MHz  is  the  source 
frequency,  £r  is  the  relative  permittivity,  jir  —  1  is  the  relative  permeability,  £0is  free  space 
permittivity,  and  [i0  is  free  space  permeability.  The  value  of  £r  —  1.15  was  chosen  to  deviate 
from  free  space  but  also  to  prevent  substantial  reflection.  Each  portion  of  the  phase  screen  is  Z(x) 
by  1  m  wide.  Incorporation  of  the  slab  into  an  FDTD  simulation  requires  use  of  an  FDTD  mesh 
cell  size  that  is  fine  enough  to  resolve  variances  in  Z(x).  Various  simulations  throughout  the 
course  of  this  research  used  resolutions  of  A/ 100  and  A/ 25,  where  A  is  the  wavelength  in  free 
space  corresponding  to  the  source  frequency  (/)  given  above.  One  consequence  of  such  a  high 
resolution  is  the  increased  FDTD  grid  size.  A  pseudo-color  rendering  of  the  dielectric  slab’s  first 
3  km  section  is  provided  in  Figure  6.  The  section  of  the  FDTD  grid  that  holds  the  slab  is  252040 
by  628  grid  cells. 
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Figure  6  0-3  km  region  of  dielectric  slab  in  red  color  from  FDTD.  The  source  wave  sweeps  from  top  to  bottom  (i.e.  flat  edge  of 

the  slab  is  exposed  first). 

3.2  Implementation  of  Near-to-Near  Field  and  Near-to-Far  Field  Transforms 

The  ionospheric  scintillation  test  case  provided  by  AFRL  for  this  project  requires  a  high- 

resolution  FDTD  grid  (Ax  =  A y  —  A/100  «  1.2  cm)  as  discussed  in  the  preceding  subsection. 

Such  a  requirement,  combined  with  the  large  size  of  the  dielectric  screen  model  (20  km)  and  the 

simulated  altitude  from  Earth’s  surface  (300  km)  resulted  in  a  number  of  difficulties  in  the 

implementation  of  the  near-to-near  and  near-to-far  field  transforms.  The  two  main  challenges 

brought  in  by  these  requirements  can  be  summarized  as: 

•  (i)  the  amount  of  increased  processing  power  needed  to  calculate  integrals  over  very 
large  contour  sizes  in  a  reasonable  time  frame,  and 

•  (ii)  the  impracticality  of  having  to  deal  with  very  large  storage  arrays  and  array  indices 
as  the  wave  propagation  time  needed  to  cover  300  km  is  over  35.7  million  FDTD  time 
steps. 

The  first  problem  was  addressed  by  taking  advantage  of  parallel  processing  capabilities  of 
modem  computer  workstations,  clusters  and  computing  clouds.  The  updates  summarized  in 
equations  (33),  (34)  and  (37)  for  the  two  transformation  techniques  were  performed  separately 
for  each  destination  point  of  interest  and  as  such,  provide  a  natural  way  of  dividing  up  the  work. 
Message  Passing  Interface  (MPI)  [27]  based  parallelization  was  used  in  this  project.  MPI 
provides  processes  that  can  be  put  to  work  in  different  portions  of  a  simulation  by  referring  to 
their  individual  identification  numbers.  These  processes  can  interact  via  special  MPI  function 
calls;  however  they  do  not  share  variables  or  arrays  in  memory.  In  the  test  case  implementation,  a 
single  MPI  process  was  used  for  FDTD  calculations,  and  63  MPI  processes  divided  up  work  for 
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the  near-to-near  and  near-to-far  field  transfonnation  calculations.  After  calculation  of  every 
FDTD  time  step,  the  single  MPI  process  broadcasted  appropriate  field  values  to  the  other  MPI 
processes.  After  broadcasting,  the  single  process  undertook  the  FDTD  calculations  for  the  next 
time  step  while  the  other  63  processes  calculated  portions  of  the  transform  with  the  received 
information.  At  the  end  of  the  simulation  the  63  processes  output  their  portions  of  the  final 
transform  data  into  a  single  data  file.  Results  obtained  from  the  parallelized  FDTD  calculation  do 
not  change  with  the  number  of  MPI  processes  used. 

The  second  challenge  was  tackled  by  utilizing  sparse  data  storage  techniques,  which  are 
elaborated  in  the  next  section.  Generally  speaking,  the  storage  technique  will  depend  on  what 
computer  language  is  used  for  the  FDTD  implementation.  Ecosystems  of  many  modem 
programming  languages  have  sparse  storage  support  available  in  one  form  or  another,  such  as 
sparse  matrices  and  arrays,  associative  arrays  and  hash  tables.  Multiple  language-agnostic 
storage  formats  such  as  Compressed  Sparse  Row,  Compressed  Sparse  Column,  and  List  of  Lists 
are  also  available. 

FORTRAN  was  the  language  of  choice  for  the  ionospheric  scintillation  test  case 
implementation.  This  language  does  not  natively  provide  an  appropriate  data  type.  The  first 
sparse  storage  fonnat  tested  with  FDTD  was  Compressed  Sparse  Row  (CSR).  The  updates 
indicated  in  equations  (33),  (34)  and  (37)  are  compatible  with  this  2-D  sparse  format  if  spatial 
coordinates  (r  or  x)  are  considered  as  one  of  the  dimensions  and  temporal  coordinates  (nn  or  t^) 
are  accepted  as  the  other  dimension.  The  CSR  storage  fonnat  divides  up  the  2-D  data  in  rows  and 
then  stores  nonzero  elements  together  with  their  column  indices.  An  auxiliary  array  indicates 
where  rows  start  and  end.  This  storage  format  is  efficient  with  various  linear  algebra  operations, 
but  random  element  retrievals  and  element  updates  require  a  linear-time  search  to  identify  the 
element  with  the  correct  column  index.  Storing  an  element  that  was  not  stored  previously 
requires  expensive  partial  array  copies.  These  issues  tend  to  degrade  performance  for  large-size 
systems. 

After  initial  tests  with  CSR  indicated  low  performance,  portions  of  the  FLIBS  project 
[28]  were  adopted  for  use  in  FDTD.  The  new  sparse  storage  support  was  thus  based  on  the  hash 
table  and  linked  list  implementations  of  the  FLIBS  project  with  various  optimizations  to  improve 
execution  speed.  An  iterator  was  also  implemented  on  the  hash  table  so  that  it  could  be  traversed 
for  diagnostics  purposes  and  loop-like  constructs  can  be  executed  with  data  from  it.  Hash  table  is 
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a  special  lookup  table  construct  that  uses  a  function  to  assign  given  key  and  value  pairs  an 
internal  location.  This  operation  can  be  performed  faster  than  anything  requiring  a  search 
regardless  of  how  large  the  data  set  becomes,  so  long  as  the  hash  function  itself  can  be  evaluated 
very  quickly  and  yields  a  reasonably  uniform  distribution  of  values  for  the  data  at  hand.  The  hash 
function  implementation  of  the  FLIBS  project  was  a  major  subject  of  the  execution  speed 
optimizations. 

3.3  Phase  Detection  Method 

Phase  detection  is  the  key  component  of  ionospheric  scintillation  simulations  in  FDTD.  The 
method  used  in  this  research  operates  by  detecting  peak  amplitude  timing  of  periodic  electric 
field  oscillations  in  FDTD  at  locations  of  interest.  The  phase  delay  of  the  incident  plane  wave  is 
calculated  by  running  the  wave  through  the  FDTD  grid.  The  maximum  values  of  Ez  for  one 
period  of  oscillation  are  found  across  the  majority  of  the  x-axis  at  a  constant  y-coordinate.  In 
other  words,  Ezmax  is  calculated  for  every  location,  ( i,j sampling )»  where  j sampling  is  located 
after  the  dielectric  slab.  The  only  i  locations  not  sampled  are  locations  close  to  the  edges  of  the 
sampling  region  due  to  the  possibility  of  edge  effects  creeping  into  the  wave  pattern.  One  period 
of  oscillation  is  defined  as  the  reciprocal  of  the  wave  frequency,  250  MHz,  which  corresponds  to 
4  ns,  or  approximately  143  time  steps  for  a  time  step  size  At  —  28  ps  (at  a  spatial  resolution  of 
Ax  —  Ay  =  A/100).  Let  the  time  step  at  which  Ezmax  occurs  at  location  (i,j  sampling)  be 
tmax (0-  The  same  procedure  is  repeated  with  an  identical  FDTD  grid  without  the  dielectric  slab. 
The  maximum  value  of  Ez  without  the  slab,  Ezmaxns,  is  calculated  across  the  x-axis  at  the  same 
y-coordinate,  j  sampling-  The  time  step  at  which  the  maximum  occurs  at  every  location, 
(i,j sampling),  to  be  identified  as,  tmaxns{i),  is  recorded  during  this  calculation.  The  phase  delay 
due  to  the  ionospheric  scintillation  process  is  calculated  as 

( Pdelay  ~  2uf  (t  maxns  (0-  tmax  (0)  (48) 

at  every  location  (i,  j  sampling)- 

4  Results  and  Discussion 

FDTD  simulation  tests  of  ionospheric  scintillation  were  performed  in  two  stages.  The  first  stage 
ensured  the  phase  screen  model  worked  as  intended.  The  second  stage  tested  the  Kirchhoff’s 
surface  integral  representation  (KSIR)  based  near-to-near  field  transform.  The  simulations  in  the 
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first  stage  used  a  higher  grid  resolution  (Ax  =  Ay  =  d/100sl.2  cm)  while  the  simulations  in 
the  second  stage  had  a  lower  grid  resolution  Ax  =  Ay  =  A/25  =  4.8  cm  to  keep  reasonable 
simulation  times  on  the  available  computing  resources. 

A  depiction  of  the  first  stage  phase  screen  testing  is  provided  in  Figure  7.  In  the  figure, 
the  incident  plane  wave  is  injected  through  the  total  field/scattered  field  (TF/SF)  contour  and 
sweeps  in  a  vertical  direction,  moving  down.  Thus  the  flat  side  of  the  phase  screen  is  illuminated 
first.  The  phase  detector  captures  electric  field  information  and  calculates  phase  delay  by 
comparing  to  a  simulation  without  the  phase  screen,  as  discussed  in  Section  3.3. 


Figure  7  Typical  configuration  of  first  stage  ionospheric  scintillation  test  with  plane-wave  source  boundary,  phase  screen  and 
phase  detector  locations  indicated.  PML  refers  to  a  perfectly  matched  layer  absorbing  medium  placed  on  the  outer  sides  of  a 
scattering  simulation  to  prevent  reflection  of  outgoing  energy. 


The  source  wave  is  a  plane  wave  with  a  sinusoidal  profile,  oscillating  at  250  MHz.  The 
wave  amplitude  is  attenuated  at  the  beginning  of  the  simulation.  It  is  ramped  up  gradually  to 
prevent  any  simulation  instabilities  due  to  transients  that  would  occur  otherwise.  The  phase 
information  is  obtained  after  the  simulation  reaches  a  steady  state.  Phase  delay  results  are  shown 
in  Figure  8  as  compared  to  the  original  data  from  AFRL.  These  results  gave  the  researchers 
confidence  to  move  forward  into  the  second  stage  of  testing. 
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Phase  Delay 


Figure  8  Near  field  phase  delay  results  comparing  phase  delay  obtained  from  FDTD  and  original  phase  delay  data  provided  by 
AFRL  for  phase  screen  generation. 

The  second  stage  tests  compared  the  near-to-near  field  transformation  results  (results 
identified  as  “test”  or  “T”  in  figures)  with  an  FDTD  simulation  that  has  two  sets  of  phase 
detectors  (results  identified  as  “compare”  or  “C”  in  figures),  adjusted  to  the  same  locations. 
While  these  simulations  are  similar  to  the  one  shown  in  Figure  7,  there  are  a  couple  key 
differences  besides  the  FDTD  grid  resolution  change.  The  first  change  is  in  the  far-field  phase 
detector  position  in  the  comparison  case.  The  near-to-near  field  transformation  contour 
(identified  as  “KSIR  contour”  in  Figure  9(a))  is  in  the  scattered  field  region  (i.e.  outside  the 
TF/SF  source  contour).  As  such,  the  far-field  phase  detector  in  Figure  9(b)  must  be  placed  in  the 
scattered  field  region  as  well.  The  second  change  is  in  the  reference  calculations  needed  for 
phase  delay  analyses.  A  scattered  field  response  needed  for  phase  detection  cannot  be  obtained 
without  a  scatterer.  As  such,  the  reference  cases  corresponding  to  simulations  in  Figure  9  include 
a  constant-height  phase  screen  instead  of  only  empty  space.  This  phase  screen  produces  a 
transmitted  field  with  uniform  phase  response. 
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Near-field  phase  detector 


Far-field  phase  detector 

J  l - j 


(a)  (b) 

Figure  9  (a)  Kirchhoff’s  surface  integral  representation  (KSIR)  based  near-to-near  field  transform  test  configuration  (results 
identified  as  “test”  or  “T”).  (b)  Corresponding  FDTD  configuration  with  two  phase  detector  arrays  (results  identified  as 
“compare”  or  “C”). 


Phase  delay  results  from  simulations  set  up  as  illustrated  in  Figure  9  are  shown  in  Figure 
10.  The  effect  of  reduced  grid  resolution  is  apparent  in  the  near-field  phase  detector  results 
(identified  as  “C  Near”).  More  important,  however  are  the  fluctuations  seen  in  the  far  field 
results.  These  suggest  that  the  phase  detection  scheme  is  not  working  well  in  the  scattered  field 
region.  Since  it  works  well  in  the  total  region,  the  next  step  taken  was  to  try  to  reconstruct  the 
total  field  from  the  sampled  scattered  field  response  by  adding  incident  field  values. 


Phase  Delay  (Ax=Ay=>./25) 


Figure  10  Phase  delay  results  from  simulations  depicted  in  Figure  9. 
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Reconstruction  of  the  total  field  requires  the  addition  of  appropriate  incident  field  values 
to  sampled  scattered  field  values.  In  a  1-D  sense,  this  reverses  the  change  done  by  equation  (24) 
of  Section  2.2  on  the  exit  side  of  the  TF/SF  source  boundary.  In  a  2-D  simulation  the  incident 
field  terms  are  extrapolated  from  a  1-D  auxiliary  simulation  for  the  incident  field  (for  details 
refer  to  [7]).  The  same  auxiliary  simulation  can  be  used  for  incident  field  reconstruction  by 
extrapolating  at  the  far-field  phase  detector  location  in  Figure  9(b),  with  results  shown  in  Figure 
1 1 .  Some  improvement  was  seen  compared  to  scattered  field-only  results.  Various  attempts  were 
made  with  heuristics  rules  aiming  to  correct  phase  delay  jumps  larger  than  2n  radians, 
smoothing  and  correction  of  the  incident  field  signal,  and  reconstruction  of  the  incident  field 
signal  as  a  sinusoidal  wave. 


Phase  Delay  (Ax=Ay=W25) 


Figure  11  Phase  delay  results  from  a  total  field  reconstruction  attempt. 


Next,  testing  efforts  were  concentrated  on  the  ability  of  using  the  phase  detection  scheme  on 
fields  sampled  in  the  total  field  region.  The  main  problem  here  was  how  to  accommodate  the 
near-to-near  field  transfonnation  into  this  picture.  Assuming  one  is  only  interested  in  the  phase  of 
the  fields  transmitted  from  the  phase  screen  and  transmissions  to  the  side  are  minimal  under  a 
small-angle  scattering  assumption,  one  can  ignore  three  sides  of  the  near-to-near  field 
transformation  contour  and  place  it  into  the  total  field  region  (Figure  12).  Edge  effects  are 
removed  as  usual  in  the  phase  delay  detection  process.  Note  that  field  amplitudes  extrapolated 
from  KSIR  are  likely  not  suitable  for  use.  These  assumptions  can  be  relaxed  with  a  different 
source  implementation,  discussed  briefly  in  the  next  section  and  illustrated  in  Figure  14. 
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TF/SF  contour 


Near-field  phase  detector 


Far-field  phase  detector 


(a)  (b) 

Figure  12  (a)  Kirchhoff’s  surface  integral  representation  (KSIR)  contour  placed  in  total  field  region  under  certain  assumptions, 
(b)  Corresponding  FDTD  configuration  with  two  phase  detector  arrays. 

The  phase  delay  results  are  shown  in  Figure  13  (a).  KSIR  produced  results  that  follow  the  phase 
changes  mostly  correctly,  except  for  a  couple  jumps.  These  can  be  readjusted  with  simple 
heuristics  (Figure  13  (b)),  since  the  phase  delay  values  cannot  be  positive  due  to  the  phase 
referencing  process  described  in  Section  3.1. 


Phase  Delay  (Ax=Ay=V25)  phase  Delay  (Ax=Ay=X/25) 


(a)  (b) 

Figure  13  (a)  Phase  delay  results  from  simulations  depicted  in  Figure  12.  (b)  After  translation  of  positive  phase  values  by  2n 
radians. 
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5  Future  Directions 

This  section  suggests  a  number  of  technical  possibilities  that  can  improve  upon  the  state  of  art 
attained  in  the  course  of  this  project.  The  ideas  herein  can  be  summarily  divided  into  two 
categories:  (i)  The  ones  that  are  meant  to  increase  numerical  accuracy  and  efficiency  of  FDTD 
simulations  for  ionospheric  scintillation  studies,  and  (ii)  the  ones  meant  to  improve  computation 
performance  and  decrease  simulation  times. 

The  discussion  in  the  previous  section  showcased  the  idea  of  using  a  near-to-near  field 
transformation  surface  within  the  total  field  region  of  a  total  field/scattered  field  (TF/SF)  source 
under  certain  conditions.  A  great  improvement  in  accuracy  can  be  attained  by  using  a  source  that 
emits  a  narrower  beam  that  lights  up  only  a  portion  of  the  ionosphere  model  in  the  simulation 
(Figure  14).  This  source  can  be  simple  or  very  sophisticated.  Choices  available  change  from 
simple  dipole  sources  [7]  to  more  sophisticated  reflectionless  dipole  sources  [13],  to  eigenmode 
injection  surfaces  [29],  to  antenna  models  with  various  degrees  of  complexity  [7].  Additionally, 
methods  such  as  “bootstrapping”,  where  the  source  wave  is  recorded  from  an  auxiliary  FDTD 
simulation  for  later  use  [7],  and  Huygens  subgridding  [30],  which  allows  combination  of  a  high- 
resolution  source  simulation  with  the  main  ionospheric  scintillation  simulation  through  virtual 
interfaces,  can  be  utilized. 


PML 


Figure  14  A  source  with  a  narrow  beam  profile  placed  into  the  KSIR  contour.  The  width  of  the  phase  screen  can  be  reduced  since 
most  of  it  is  not  illuminated  by  the  source. 
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It  is  possible  to  have  broadband  phase  detennination  if  it  is  perfonned  in  the  frequency 
domain  via  Discrete  Fourier  Transfonn  (DFT).  More  details  on  efficient  DFT  implementation  in 
FDTD  can  be  found  elsewhere  [13].  DFT  transformed  results  will  be  complex  numbers  with 
amplitude  and  phase  for  every  frequency  covered.  Since  DFT  implies  the  electric  field  signal 
sampled  from  FDTD  to  be  periodic  or  quasi-periodic  in  nature,  this  method  requires  some 
changes  to  the  source  wave  (i.e.  any  non-periodic  source  wave  patterns  have  to  be  repeated 
multiple  times  during  the  course  of  the  simulation).  Additionally,  windowing  and  other 
considerations  with  Fourier  transforms  might  apply. 

Beach  et  al.  [1]  performed  a  sparse  sampling  study  of  the  amplitude  scintillation  index 
(S4).  This  study  provides  important  criteria  needed  to  get  a  good  approximation  of  S4  in  the  event 
of  inadequate  sampling.  In  this  vein,  practical  benefit  might  arise  if  FDTD  simulation  results  are 
under-sampled  prior  to  S4  calculations.  Filtering,  smoothing  or  other  sophisticated  signal 
correction  heuristics  can  also  be  applied.  Such  techniques  are  commonly  used  with  FDTD  for 
certain  cases,  for  example  a  Pade  approximant  based  spectral  analysis  procedure  can  be  used  for 
acquisition  of  sharp  resonances  [7]. 

Perfonnance  of  the  simulations  can  be  improved  by  using  graphics  accelerators.  Graphics 
accelerators,  originally  developed  for  sophisticated  scene  rendering,  work  well  on  data  sets  with 
regular  strides  and  access  patterns.  Many  articles  on  FDTD  acceleration  reporting  excellent 
results  can  be  found  in  literature  [7].  The  ionospheric  scintillation  simulation  can  be  accelerated 
similarly  if  the  hash  table/linked  list  design  is  replaced  with  a  more  cooperative  one.  Since  the 
near-to-near  and  near-to-far  field  transform  implementations  have  predictable  data  access 
patterns,  it  should  be  possible  to  design  and  implement  a  different  data  structure  that  is  more 
compatible  with  graphics  accelerators. 

Last  but  not  least,  a  recently  developed  FDTD  variant  [29,31,]  allows  calculation  of 
stochastic  mean  and  variance  estimates  with  FDTD.  If  this  method  can  be  extended  to  address 
the  problem  at  hand,  stochastic  modeling  of  ionospheric  scintillation  would  provide 
unprecedented  time  savings  by  generating  statistical  information  that  can  be  obtained  from 
hundreds  or  thousands  of  randomized  FDTD  models.  Such  statistical  infonnation  is  needed  to 
assess  or  counter  effects  of  scintillation  on  communications,  such  as  channel  fading. 
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6  Conclusion 


This  study  investigated  FDTD  modeling  of  ionospheric  scintillation  effects.  Methods  of  near-to- 
near  field  and  near-to-far  field  transformation  were  implemented  to  enable  FDTD  simulation  of 
scintillation  perturbed  signals  at  large  distances.  These  methods  were  augmented  by  hash  table 
based  sparse  data  storage  and  parallelized  evaluation  to  make  them  usable  for  an  AFRL  provided 
test  case  scenario. 

Research  undertaken  evaluated  phase  detection  in  FDTD  and  its  compatibility  with  near- 
to-near  field  and  near-to-far  field  transformation  techniques  needed  for  extrapolation  based  range 
extension.  The  results  obtained  indicate  it  is  possible  to  use  these  techniques  together  with  a 
different  source  implementation.  Finally,  a  discussion  of  future  research  directions  and  potential 
improvements  to  the  current  study  was  included.  FDTD  has  excellent  potential  in  terms  of 
simulation  accuracy  improvements  and  performance  improvements  on  modern  computing 
hardware.  This  is  besides  its  potential  for  stochastic  modeling,  possibilities  for  advanced  antenna 
and  ionospheric  scintillation  models  and  removal  of  various  restrictions  found  in  the  phase 
screen  method.  It  is  the  researchers’  opinion  that  future  research  on  FDTD-based  ionospheric 
scintillation  modeling  will  bring  exciting  results  far  beyond  the  reach  of  the  current  state  of  art. 
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